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STOKES EQUATIONS 
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Abstract. This paper introduces a weak Galerkin (WG) finite element method for the Stokes 
equations in the primary velocity-pressure formulation. This WG method is equipped with stable 
finite elements consisting of usual polynomials of degree fc > 1 for the velocity and polynomials 
of degree k — 1 for the pressure, both are discontinuous. The velocity element is enhanced by 
polynomials of degree fc — 1 on the interface of the finite element partition. All the finite element 
functions are discontinuous for which the usual gradient and divergence operators are implemented 
as distributions in properly-defined spaces. Optimal-order error estimates are established for the 
corresponding numerical approximation in various norms. It must be emphasized that the WG finite 
element method is designed on finite element partitions consisting of arbitrary shape of polygons or 
polyhedra which are shape regular. 
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1. Introduction. In this paper, we are concerned with the development of weak 
Galerkin (WG) finite element methods for the Stokes problem which seeks unknown 
functions u and p satisfying 

(1.1) -Au + Vp = f inn, 

(1.2) V ■ u = inn, 

(1.3) u = g on dn, 

where is a polygonal or polyhedral domain in M'' (d = 2, 3). 

The weak form in the primary velocity-pressure formulation for the Stokes prob- 
lem ([rT])~([r3l) seeks u e [H^{n)]'^ and p £ Ll{n) satisfying u = g on dn and 

(1.4) (Vu,Vv)-(V-v,p) = (f,v), 

(1.5) (V-u,g)=0, 

for all V G [i?Q(r2)]'' and q G Lo(^)- The conforming finite element method for (jl.ip - 
(jl.3|) developed over the last several decades is based on the weak formulation (|1.4p - 
(jl.Sp by constructing a pair of finite element spaces satisfying the inf-sup condition 
of Babuska [1] and Brezzi [3]. Readers are referred to [6] for specific examples and 
details in the classical finite element methods for the Stokes equations. 

Weak Galerkin refers to a general finite element technique for partial differential 
equations in which differential operators are approximated by weak forms as distri- 
butions for generalized functions. Thus, two of the key features in weak Galerkin 
methods are (1) the approximating functions are discontinuous, and (2) the usual 
derivatives are taken as distributions or approximations of distributions. A weak 
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Galerkin method was first introduced and analyzed for second order elliptic equations 
in and in a conference in Nankai University in the summer of 2011 by one of the 
authors. The objective of this paper is to develop a weak Galerkin finite element 
method for (|l.ip - (|1.3[) that is efficient and robust by allowing the use of discontinuous 
approximating functions on finite element partitions consisting of arbitrary polygons 
or polyhedra with certain shape regularity. 

In general, weak Galerkin finite element formulations for partial differential equa- 
tions can be derived naturally by replacing usual derivatives by weakly-defined deriva- 
tives in the corresponding variational forms, with the option of adding a stabilization 
term to enforce a weak continuity of the approximating functions. For the Stokes 
problem (jl.l|) - (ll.3p interpreted by the variational formulation (|1.4p - (|1.5p . the two 
principle differential operators are the gradient and the divergence operator defined 
in the Sobolev space [H^{n)]'^. Formally, our weak Galerkin method for the Stokes 
problem would take the following form: Find u/i and ph from properly-defined finite 
element spaces satisfying 



for all test functions v and q in test spaces. Here V^, is a weak gradient and V^- 
is a weak divergence operator to be detailed in this study (see Section [3] for details) . 
The bilinear form s(-,-) in (|1.6p is a parameter free stabilizer that shall enforce a 
certain weak continuity for the underlying approximating functions. The use of totally 
discontinuous functions and weak derivatives in the WG formulation provides the 
numerical scheme with many nice features. First, the construction of stable elements 
for the Stokes equations under WG formulation is straight forward with standard 
polynomials. Secondly, the WG method allows the use of finite element partitions with 
arbitrary shape of polygons in 2D or polyhedra in 3D with certain shape regularity. 
The later property provides a convenient and useful flexibility in both numerical 
approximation and mesh generation. Thirdly, our WG formulation is parameter- 
free and has competitive number of unknowns since lower degree of polynomials are 
used on element boundaries, and the unknowns corresponding to the interior of each 
element can be eliminated from the system. 

The paper is organized as follows. In Section [21 we introduce some standard 
notations in Sobolev spaces. Two weakly-defined differential operators, weak gradient 
and weak divergence, arc introduced in Section [H The WG finite element scheme 
for the Stokes problem (|l.ip - (|1.2p is developed in Section |4l In Section [5l we shall 
study the stability and solvability of the WG scheme. In particular, the usual inf- 
sup condition is established for the WG scheme. In Section \6\ we shall derive an 
error equation for the WG approximations. Optimal-order error estimates for the 
WG finite element approximations are derived in Section [7] in virtually an norm 
for the velocity, and norm for both the velocity and the pressure. In Section [8l 
we make some concluding remarks by mentioning some outstanding issues for future 
consideration. Finally, we present some technical estimates for quantities related to 
the local projections into various finite element spaces in Appendix [Al 

2. Preliminaries and Notations. Let D be any open bounded domain with 
Lipschitz continuous boundary in M*^, d = 2, 3. We use the standard definition for the 
Sobolev space H^{D) and the associated inner product (•,-)s,_d, norm || • ||s,_d, and 
seminorm | • \s^d for any s > 0. For example, for any integer s > 0, the seminorm 



(1.6) 
(1.7) 



■) - (V^ • V,ph) + s(Uft,v) = (f,v), 
(Vu, ■VLh,q) = Q 
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I • U_D is given by 
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with the usual notation 



d 



a = (ai, . . . ,ad), 



a\ 



The Sobolev norm \\ ■ \\m,D is given by 



EH 



3=0 
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The space H^{D) coincides with L^{D), for which the norm and the inner product 
are denoted by \\-\\d and {■,-)d, respectively. When D = 57, we shall drop the subscript 
D in the norm and inner product notation. 

The space H{div; D) is defined as the set of vector-valued functions on D which, 
together with their divergence, are square integrable; i.e.. 



3. Weak Differential Operators and Their Approximations. The key to 

weak Galerkin methods is the use of weak derivatives in the place of strong derivatives 
that define the weak formulation for the underlying partial differential equations. The 
two differential operators used in the weak formulation (|1.4p and (|1.5p are gradient 
and divergence. Thus, it is essential to introduce a weak version for both the gradient 
and the divergence operator. In [10], a weak divergence operator has been introduced 
and employed to the mixed formulation of second order elliptic equations. In |9] 
and 0, a weak gradient operator was introduced for scalar functions. Those weakly 
defined differential operators shall be employed to the Stokes problem ()1.4|) - (|1.5p in 
a weak Galerkin approximation. For convenience, the rest of the section will review 
the definition for the weak gradient and the weak divergence, respectively. Note that 
the weak gradient shall be applied to each component when the underlying function 
is vector- valued, as is the case for the Stokes problem. 

3.1. Weak gradient and discrete weak gradient. Let K be any polygonal 
or polyhedral domain with boundary dK. A weak vector-valued function on the 
region K refers to a vector- valued function v — {vq, v;,} such that Vq G [L'^{K)Y and 
Vf, G [H^{dK)Y. The first component Vq can be understood as the value of v in K, 
and the second component Vf, represents v on the boundary of K. Note that Vf, may 
not necessarily be related to the trace of Vq on dK should a trace be well-defined. 
Denote by V{K) the space of weak functions on K] i.e., 



i/(div; D)^{w. V e [L'{D)Y, V • v G L'iD)] . 



The norm in iJ(div; D) is defined by 




(3.1) V{K) = {v = {vo, Vb} : vq G [L\K)Y, G [i/^ {dK)Y}. 
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The weak gradient operator is defined as follows. 

Definition 3.1. The dual of [Li^i^Kyf can he identified with itself by using the 
standard inner product as the action of linear functionals. With a similar inter- 
pretation, for any v € V{K), the weak gradient ofv is defined as a linear functional 
VujV in the dual space of[H{div,K)Y whose action on each q € [H{div,K)]''- is given 
by 

(3.2) (V^v, q)K = -(vo, V ■q)K + (vft, q ■ n)dK, 

where n is the outward normal direction to dK, (vq, V • q)K = Sk^^^ ' l)^^ 
the action of vq on V • q, and (v;,, q ■ n)gK = Jqj^ vtQ ■ nrfs is the action of q ■ n on 

The Sobolev space [H^ (K)]''' can be embedded into the space V{K) by an inclusion 
map iv ■■ [H^{K)]'^ ^ V{K) defined as follows 

iv(0) = {'/'k,<^|o^}, <Pe[H\K)f. 

With the help of the inclusion map ly, the Sobolev space [H^{K)]'^ can be viewed as 
a subspace of V{K) by identifying each </> e [H^{K)]'^ with i\>{(j)). 

Let Pr{K) be the set of polynomials on K with degree no more than r. 

Definition 3.2. The discrete weak gradient operator, denoted by ^w.r,K, is 
defined as the unique polynomial (Vw,r,i<:v) € [Pr(isr)]''^'^ satisfying the following 
equation, 

(3.3) {Vy,,r,KV,q)K = -{^0,^ ■q)K + {vb,q-n)dK, Vq e [Pr{K)f'"^. 



3.2. Weak divergence and discrete weak divergence. To define weak di- 
vergence, we require weak function v = {vcVft} such that vq G [L^(ii')]'^ and 
Vf, • n e H~^{dK). Denote by V{K) the space of weak vector-valued functions 
on K] i.e., 

(3.4) V{K) = {v = {vo, Vft} : vq S [L'{K)]\ • n G H-"^{dK)}. 

A weak divergence operator can be defined as follows. 

Definition 3.3. The dual of L'^{K) can be identified with itself by using the 
standard inner product as the action of linear functionals. With a similar inter- 
pretation, for any v G V{K), the weak divergence ofv is defined as a linear functional 
Vw ■ V in the dual space of {K) whose action on each G {K) is given by 

(3.5) (V,„ • V, ip)K = -(vo, V<^)/f + (V6 ■ n, 'p)dK, 

where n is the outward normal direction to dK, (vq, V(/j)i<- is the action ofvo on V</5, 
and (vb • n,ip)QK is the action of ■ n on G H~{dK). 

The Sobolev space [H^{K)Y can be embedded into the space V{K) by an inclusion 
map iv : [H\K)Y ^ V{K) defined as follows 

iv{4>) = {4>\K, <k\dK], <k(^[H\K)\''. 
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Definition 3.4. A discrete weak divergence operator, denoted by ^w.r,K'j 
defined as the unique polynomial {^w.r,K ' v) G Pr{K) that satisfies the following 
equation 

(3.6) C^w.r^K ■v,'p)k = -{va,'^^)K + {vb ■'n,ip)dK, \f(p e Pr{K). 

4. A Weak Galerkin Finite Element Scheme. Let Th be a partition of the 
domain n with mesh size h that consists of arbitrary polygons/poly hedra. Assume 
that the partition Th is WG shape regular - defined by a set of conditions as detailed 
in [TU] and [5]. Denote by £h the set of all edges/flat faces in Th, and let £^ = £h\dil. 
be the set of all interior edges/faces. 

For any integer fc > 1, we define a weak Galerkin finite element space for the 
velocity variable as follows 

Vh = {w^ {vo,v,} : {vo,v,}|t G [Pk{T)f x [Pfc-i(e)]^ e C ^T} . 

We would like to emphasize that there is only a single value v;, defined on each edge 
e ^ £h- For the pressure variable, we have the following finite element space 

Wh^{q: qeLlin), gjr G Pfc-i(T)} . 

Denote by the subspace of Vh consisting of discrete weak functions with vanishing 
boundary value; i.e., 

= {v = {vo, vj G Vh, Vf, = on dn} . 

The discrete weak gradient Vu),fc-i and the discrete weak divergence {'Vw,k-i') on the 
finite element space Vh can be computed by using (|3.3p and (|3.6p on each element T, 
respectively. More precisely, they are given by 

(V„,fc_iv)|T = V^,fc_i,T(v|T), V V G Vh, 

(Vn,,fc_i • v)|t = Vu,,fc-i,T • (v|t), V V G Vh. 

For simplicity of notation, from now on we shall drop the subscript fc — 1 in the 
notation ^w,k-i and (S/w,k-i') for the discrete weak gradient and the discrete weak 
divergence. The usual inner product can be written locally on each element as 
follows 

(Vt„v, V.ujw) = ^ (V^„v, V^w)t, 

Ten 

(V„ • V, (?) = ^ (V,„ • V, q)T. 
Ten 

Denote by Qo the projection operator from [L^(r)]'^ onto [Pk{T)]'^. For each 
edge/face e G £h, denote by Qb the projection from [L^{e)]'^ onto [Pk-i{e)]'^. We 
shall combine Qq with Qi, by writing Q}^ = {Qq, Qb}- 

We are now in a position to describe a weak Galerkin finite element scheme for 
the Stokes equations (ll.ll) - (|1.3p . To this end, we first introduce three bilinear forms 
as follows 

s(v, w) = ^ h'j}{QbWo-^fb, Qb^Q-^b)QT^ 
Ten 

a(v, w) = (V^,v, V^,w) + s(v, v^f), 
6(v, q) = (Vu, • V, q). 
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Weak Galerkin Algorithm 1. A numerical approximation for U.1\) - U.S\) 
can be obtained by seeking U/, = {uo,Ut} G Vh and £ Wh such that Uf, — QbS on 
dfl and 

(4.1) a{uh, v) - 6(v, ph) = if, vo), 

(4.2) b{uh,q)^0, 

for all V = {vq, Vb} G Vf^ and q G Wh- 

5. Stability and Solvability. The WG finite element scheme (|4TT|) - (|4?2l) is a 
typical saddle-point problem which can be analyzed by using the well known theory 
developed by Babuska [T] and Brezzi [5]. The core of the theory is to verify two 
properties: (1) boundedness and a certain coercivity for the bilinear form a(-, •), and 
(2) boundedness and inf-sup condition for the bilinear form •). 

The finite element space is a normed linear space with a triple-bar norm given 

by 

(5.1) lllvf = Yl + E '^r'WQb^o - MIt- 

Ten Ten 

We claim that ||| • ||| indeed provides a norm in V^. For simplicity, we shall only verify 
the positive length property for ||| • |||. Assume that |||v||| = for some v G V,". It 
follows that 

= (V„v, V^v) -I- ^ h'j} {QbVn - ^^b, Qb^o ~^fb)QT^ 
TeTh 

which implies that Vu,v = on each element T and QbVo = Vb on dT. Thus, we have 
from the definition that for any r G [Pk-i{T)Y'^'^ 

0= (V^,v,t)t 
= -(vo, V • t)t + {vb^T ■ n)dT 
= (Vvo, t)t - (vo - Vb, r • n)dT 
= (Vvo, r)T - (QbVo - Vb, T • n)gT 
= (Vvo,r)T. 

Letting r = Vvq in the equation above yields Vvq = on T G Th- It follows that 
Vg = const on every T € Th- This, together with the fact that QbVo = Vb on dT and 
Vb = on implies that Vq = and Vb = 0. 

Note that ||| • ||| defines only a semi-norm in Vh- It is not hard to see that a(v, v) = 
|||v||| for any v G Vh. In fact, the trip-bar norm is equivalent to the standard H^- 
norm, but was defined for weak finite element functions. It follows from the definition 
of III • III and the usual Cauchy-Schwarz inequality that the following boundedness and 
coercivity hold true for the bilinear form a(-, •). 

Lemma 5.1. For any v,w G V^, we have 

(5.2) |a(v,w)| < |||v||||||w|||, 

(5.3) a(v,v) = lllvf . 
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In addition to the projection — {Qq, Qb} defined in the previous section, let 
Qh and Q/i be two local projections onto Pk-i{T) and [Pk-i{T)Y^'^ , respectively. 

Lemma 5.2. The projection operators Qu, Q^, and Qh satisfy the following 
commutative properties 

(5.4) V„(g,,v) = Q^Vv), Vve[H\n)f, 

(5.5) V^-(0;,v)=Q,,(V-v), V vei/(div,fi). 

Proof. Using p.3p . we have 

C^w{Qh-v), q)T = -(Qov, V • g)T + (QbV, q ■ n)oT 

for all q G [Pfe_i(r)]'*^'*. Next, we use the definition of Qh and Qh and the usual 
integration by parts to obtain 

-(Qov, V • (7)t + (QbV, q • n)oT = -(v, V ■ (7)t + (v, q ■ n)gT 

= (Vv, q) 
= (Q,.(Vv), q). 

Thus, 

(v^(g„v), g)T = (Q/.(Vv), g), vqe [Pfc_i(r)]'^x^ 

which verifies the identity (|5.4p . 

To verify (|5.5I) . we use the discrete weak divergence p.6p to obtain 

(V^ • (Q/iv), (p)t = -(Qov, V(p)t + (QbV • n, ip)dT 

for all (/3 G Pk-i{T). Next, we use the definition of Qh and Q/i and the usual integra- 
tion by parts to arrive at 

-(Qov, V<^)t + (QbV • n, Lp)oT = -(v, V<^)t + (v • n, Lp)oT 

= (V • V, (p)t 
= (Qh(V-v), ^)r. 

It follows that 

(Vu, • (Q,,v), ^)t = (Q/,(V • v), (^)t, V G Pfe-i(T). 

This completes the proof of (|5.5p . and hence the lemma. □ 

For the bilinear form b{-, ■), we have the following result on the inf-sup condition. 
Lemma 5.3. There exists a positive constant (3 independent of h such that 

(5.6) sup > (3\\p\\ 
for all p G Wh ■ 

Proof. For any given p G Wh C Ll{^l), it is well known [21 HI [Sj El [7] that there 
exists a vector- valued function v G [i/Q(r2)]'^ such that 
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where C > is a constant depending only on the domain Q. By setting v — QhV G V/i, 
we claim that the following holds true 

(5-8) |||v||| < Co||v||i 

for some constant Co- To this end, we use equation (|5.4p to obtain 

(5.9) ii^-^iIt = E = iiQ/'Vviii < iivvf . 

TeTh TGTh TeTh 

Next, we use (|A.4p . (lA.ll) . and the definition of Qb to obtain 

= E hT'WQbiQo^-mlr 
TeTh 

< E /^t'IIQov-vIIIj, 
Ten 

<CJ2 {hr'WQo^ - Ht + l|V(Qov - v)|||) 

Ten 

(5.10) <C||Vvf. 

Combining the estimate (|5.9p with (|5.10p yields the desired inequality 
It follows from (|5.5I) and the definition of Qh that 



6(v, p) = (V„-(Q,^v), p) = (Q^V • v), p) = (V • V, p). 
Using the above equation, (15. 7p and (|5.8p . we have 



ll|v||| C'o||v||i 

for a positive constant /3. This completes the proof of the lemma. □ 

It follows from Lemma 15.11 and Lemma 15.31 that the following solvability holds 
true for the weak Galerkin finite element scheme (j4.ip - (14.21) . 

Lemma 5.4. The weak Galerkin finite element scheme ^.l^ - ^T^ has one and 
only one solution. 

6. Error Equations. Let u/i = {uo,Ub} G Vh and ph G Wh be the weak 
Galerkin finite element solution arising from the numerical scheme (|4.ip - (|4.2p . De- 
note by u and p the exact solution of (|l.ip - (|1.3l) . The projection of u in the finite 
element space Vh is given by 

QhU = {Qou, Qfcu}. 

Similarly, the pressure p is projected into Wh as QhP- Denote by e/i and eh the 
corresponding error given by 

(6.1) eh = {eo, et} ^ {QqU-Uq, QbU-Ub}, Sh = QhP - Ph- 
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The goal of this section is to derive two equations for which the error e/j and Sh shall 
satisfy. The resulting equations are called error equations, which play a critical role 
in the convergence analysis for the weak Galerkin finite element method. 

Lemma 6.1. Let (w;p) G [H^{il)]'^ x L^(f2) be sufficiently smooth and satisfy the 
following equation 

(6.2) -Aw + Vp = ?7 

in the domain il. Let QhVf — {Qqw, Qbw} and QhP be the L^ projection of (w; p) 
into the finite element space Vh x Wh ■ Then, the following equation holds true 

(6.3) (V^(gftw), V„v) - (V^-v,Q,,p) = (?7,vo) +£w(v) - 0p(v) 
for all V e V;^ , where £w o,nd Op are two linear functionals on V/^ defined by 

^p(v) = ^ {^0- Vb, (p - Qhp)n)9T- 

Proof. First, it follows from (|5.4p . p.3p . and the integration by parts that 

(V.^(gftw),V„v)T = (Q?.(Vw),V^v)t 

= -(vo, V • Q/j(Vw))t + (vt, Q,i(Vw) ■ n)oT 
= (Vvo, Q/i(Vw))t - (vo - V6, Q/i(Vw) ■ n)aT 

(6.4) = (Vw, Vvo)t " (vo - Vb, Q,,(Vw) • n)^^,. 

Next, by using (|5.5p . p.6p . the fact that J^TeTh^^^'P n)^^ = and the integration 
by parts, we obtain 



(V^-v,Q^p) (vo, V(Q,,p))t + J2 {^b,iQhp)n} 



dT 



^ (V • vo,Q/,./9)t - ^{'^Q- vf,, (Q/,.p)n)aT 

^ (V • Vo, p)t - ^ (vo - Vf,, {<^hP)^)dT 

- X! ('^o, Vp)t + ^ (vo,pn)aT - ^ i^o - v^, (Q,ip)n)aT 
TeTh TeTh TeTh 

Ten TeTh Ten 

-(vo, Vp) + X! '^^0 ^ ^ Q/jP)n)aT, 



which leads to 



(6.5) (vo, Vp) = -(V^-v, Qhp) + ^{^0- vb, (p - Qhp)n)aT- 

Ten 

Next we test (|6.2I) by using Vo in v = {vo, V;,} £ to obtain 

(6.6) - (V ■ (Vw), Vo) + (Vp, vo) = iv, vo). 
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It follows from the usual integration by parts that 

-(V • (Vw), vo) = ^ (Vw, Vvo)t - ^ (vo-vfc, Vw • n)^^, 
TeTh Ten 

where we have used the fact that X^TeTh^^f"^^ ■ ^)dT — 0. Using (|6.4p and the 
equation above, we have 

- (V • (Vw), vo) = (V„(Q,,w), V.^v) 

(6.7) - ^ {'^0- Vf,, Vw ■ n - Q/,.(Vw) • n)^^. 

Ten 

Substituting and (HJ]) into dM]) yields 

(V„(Q,,w), V^v) - (V„-v, Qhp) = iv, Vo) + C(v) - 9p{v), 

which completes the proof of the lemma. □ 

The following is a result on the error equation for the weak Galerkin finite element 
scheme (|4A|) - (|42l) . 

Lemma 6.2. Let and Sh be the error of the weak Galerkin finite element 
solution arising from ^.lj^ - (J7B^ , as defined by id. Then, we have 

(6.8) a{eh, v) - 6(v, en) = Vu,j9(v), 

(6.9) b{eh, q) = 0, 

for all V S and q G Wh, where </?u,p(v) = ^u(v) — 0p(v) + s{QhU,v) is a linear 
functional defined on . 

Proof. Since (u;p) satisfies the equation (|6.2p with r/ ~ i, then from Lemma [01 
we have 

(V^(g,,u),V^v) -(V^-v,Q^p) = (f,vo)+^u(v) -0p(v). 
Adding s{QhU,v) to both side of the above equation gives 

(6.10) a{QhU, v) - 6(v, QhP) = (f , vq) + ^u(v) - 9p{v) + s(g,,u, v). 
The difference of (|6.10p and (|4.ip yields the following equation, 

a{eh, v) - 6(v, Sh) = ■^u(v) - 9p{v) + s{QhU, v) 

for all V e V,°, where = {gq, ef,} = {QqU - Uo, QbU - Uf,} and Eh = QhP - Ph- 
This completes the derivation of ()6.8p . 

As to (16.91) . we test equation (|1.2I) by q e and use (j5.5p to obtain 

(6.11) 0= (V-u,g) = (V„ •Q,,u,(7). 
The difference of (|6.1ip and (|4.2p yields the following equation 

6(e/i,g) = 

for all q £ VF/i. This completes the derivation of (16. 9p . □ 
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7. Error Estimates. In this section, we shall establish optimal order error es- 
timates for the velocity approximation in a norm that is equivalent to the usual 
iJ^-norm, and for the pressure approximation ph in the standard norm. In addi- 
tion, we shall derive an error estimate for in the standard norm by applying 
the usual duality argument in finite element error analysis. 

Theorem 7.1. Let (u;p) e [Hl{Vt) r\H^+''[nyf x {Ll{n) nH''{n)) with k > 1 
and {uh;ph) G Vh x Wh be the solution of i f J. j|) -/ f77igj) and respectively. 
Then, the following error estimate holds true 

(7.1) |||g,.u - u,,||| + WQhP - PhW < Ch\\\u\\k+i + \\p\\k). 



Proof. By letting v ~ eii in (|6.8p and q — Eh in (|6.9p and adding the two resulting 
equations, we have 

(7-2) ||e,jf = (pu,p{eh)- 

It then follows from (IA.6I) - (|A.8|) (see Appendix A) that 

(7.3) |||e,f <C/i'=(||u||fe+i + |HU)|||e^|||, 

which implies the first part of (|7.ip . To estimate ||£;i||, we have from (|6.8I) that 

b{v,eh) = a(e,i,v) - <y5u,p(v). 
Using the equation above, (15. 2p . (|7.3p and (|A.6|) - (|A.8p . we arrive at 

|6(v,£^)| <c;.^(||u||fe+i + |biu)|||v|||. 

Combining the above estimate with the inf-sup condition (j5.6p gives 

\\eh\\<ChH\Mk+i + \\p\\k), 
which yields the desired estimate (|7.1I) . □ 

In the rest of this section, we shall derive an L^-error estimate for the velocity 
approximation through a duality argument. To this end, consider the problem of 
seeking (■(/'; such that 

(7.4) - Ai/; -I- Ve = eo in n, 

(7.5) V-V' =0 inf], 

(7.6) V = ondVL. 

Assume that the dual problem has the [H'^{^l)Y x i/^(il)-regularity property in the 
sense that the solution (t/";^) G [H'^{VI)Y x H^{^) and the following a priori estimate 
holds true: 

(7.7) 11^112 + ||eili<C||eo||. 

Theorem 7.2. Let (u;p) e [H^{n) n H'^+^Q)]'^ x {Ll{n) n H''{n)) with k > 1 
and (uii'jPh) G x Wh be the solution of j|) - fO)) and respectively. 
Then, the following optimal order error estimate holds true 

(7.8) IIQou - uoll < C7/i'=+i(||u|U+i + \\p\\k). 
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Proof. Since (V'iO satisfies the equation (16. 2p with = Gq = QqU — Uq, then from 
(|6.3p we have 

(V^Q^V, V^v) - (V,„ ■ V, QhO = (eo, vo) + ^^.(v) - %(v), Vv e 
In particular, by letting v = we obtain 

lleolP = (Vu,Q,.V', V^e,,) - (V^ • eh,QhO - ^v(eh) + %(e/.). 
Adding and subtracting s{Qfi4',ef,) in the equation above yields 



where <^^,5(v) = ^^(e,,) - 6*^(6,1) + s{Q , et) . It follows from ((O)) . (173)) and ((OTI) 
that 

Combining the above two equations gives 

||eo|p ^ a{eh,Qhip) - KQh-^P, £h) - 'P^A^h)- 
Using (j6.8p and the equation above, we have 
(7-9) ||eo||^ = Vu.piQhi') - V4:.d^h)- 



To estimate the two terms on the right hand side of (|7.9p . we use the inequalities 
(IA.6I) - (|A.8|) with (w; p) = v = e/i, and r = 1 to obtain 



(7.10) 



\v^,deh)\ < ChiUh + m\i)\\\eh\\\ < Ch\\\e4\ lleoll. 



where we have used the regularity assumption (|7.7p . Each of the terms in (p-a,piQh'>P) 
can be handled as follows. 

(i) For the stability term s{QfiU, Qh'ip), we use the definition of Qi, and (jA.4[) to 
obtain 



\s{QhU, Qh^)\ 



hj,^{Qb{Qou- u), Qb{Qo-ip - ■4'))aT 



1/2 



<C/l'=+^||u|U+i||V||2. 



(ii) For the term iuiQh'ip), we first use the definition of Qt and the fact that 
'0 = on dft to obtain 

E (V' - Qbi^, Vu • n - Q/i(Vu) • n)^^, = ^ 0^ - Qbi^' ' "■)aT = 0. 
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Thus, 



J2 (Qoi' - Qb^, Vu • n - Q,,(Vu) • n) 



OT 



Ten 



1/2 



IQoV'-V'll 



dT 



<ch'+'\\u\\k+i\m,. 

(iii) For the term 9p{Qfiip), we first use the definition of Qi, and the fact that 
■0 = on dfl to obtain 

E (V' - Qbi', {p - QhP)n)aT = Y Qbi^,pn)dT = 0. 
Ten Ten 

Thus, from (IA.4I) and (jA.Sp we obtain 



9piQh4')\ 



Ten 



Ten 



< 



E hT\\p-QhP\\lT] f E hr'WQoi'-HlT] 
''en I \Ten I 



1/2 



\Ten 



<c/^''^+^lb|U.||V^||2. 

The three estimates in (i), (ii), (iii), and the regularity (|7.7p collectively yield 



(7.11) 



\^^AQh^)\ < ch'^+WMk+i + \\p\\k)\m2 

<C/i^+i(||u|U+i + |HU)||eo||. 



Finally, substituting ((7?TU| and ([7111) into ([7^ gives 

lleof < Ch'^+WMk+i + |b||fe)||eo|| +C/i|||e„,||| ||eo||. 

It follows that 

\\e4<Ch''+\\\uU+i + \\p\\k) + Ch\\\en\l 
which, together with Theorem 17. 1[ completes the proof of the theorem. □ 



8. Concluding Remarks. This paper introduced a new finite element method 
for the Stokes equations by using the general concept of weak Galerkin. The scheme 
is applicable to finite element partitions of arbitrary polygon or polyhedra. The paper 
has laid a solid theoretical foundation for the stability and convergence of the weak 
Galerkin method. There are, however, many open issues that need to be investigated 
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in future work. Here we would like to list a few for interested readers to consider: (1) 
how the discretized linear systems can be solved efficiently by using techniques such 
as domain decomposition and multigrids? (2) can the weak Galerkin scheme for the 
Stokes equations be hybridized? If so, how such a hybridization may help in variable 
reduction and solution solving? and (3) what superconvergence can one develop for 
the weak Galerkin method? (4) is the weak Galerkin method more competitive than 
other existing finite element schemes in practical computation? (5) what stability do 
weak Galerkin methods have in other norms such as LP,p > 1? 

Appendix A. 

In this Appendix, we shall provide some technical results regarding approximation 
properties for the projection operators Qh, Q/i, and Qh- These estimates have been 
employed in previous sections to yield various error estimates for the weak Galerkin 
finite element solution of the Stokes problem arising from the scheme (j4.1[) - (14.21) . 

Lemma A.l. Let Th be a finite element partition of D, satisfying the shape reg- 
ularity assumption as specified in flOf and w E [H^^^(n)]'^ and p G H^'{n) with 
1 < r < k. Then, for < s < 1 we have 

(A.l) '^T l|w - Qo^Wh < h'^^+'^Ml^,, 

Ten 

(A.2) J2 ''tW^^ - Qh{Vw)\\l, < Ch'-\\w\\l+„ 

TeTh 

(A.3) hl^\\p-Q,p\\l,<Ch^^\\p\\l. 
Ten 

Here C denotes a generic constant independent of the meshsize h and the functions 
in the estimates. 

A proof of the lemma can be found in [lOj . which is based on some technical 
inequalities for functions defined on polygon/polyhedral elements with shape regular- 
ity. We emphasize that the approximation error estimates in Lemma lA.ll hold true 
when the underlying mesh Th consists of arbitrary polygons or polyhedra with shape 
regularity as detailed in [10] and [8]. 

Let T be an element with e as an edge/face. For any function g G H^{T), the 
following trace inequality has been proved to be valid for general meshes satisfying 
the shape regular assumptions detailed in [TO] : 

(A.4) \\9\\l<C{h^%rT + hT\\S7gr^). 

Lemma A. 2. For any v — {vo, vt,} G Vh, we have 

(A.5) l|Vvo|||^<C|||t;f. 

Ten 

Proof. For any v = {vo,vt,} G Vh, it follows from the integration by parts and 
the definitions of weak gradient and Qb, 

(Vvo, Vvo)t = -(vo, V • Vvo)t + (vq, Vvo ■ n)^,, 

= -(vo, V • Vvo)t + (vfc, Vvo • n)gj, + (vq - v,,, Vvq • n)gj, 
= (Vt„v, Vvo)t + (QbVo - Vb, Vvq • n)gj,. 
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By applying the trace inequality (|A.4p and the inverse inequality to the equation 
above, we obtain 



iVvolll < C(||V^v||t||Vvo||t + h^' HQbVo - V6||aT|| Vvo||t) 



Thus, 



iVvolll < C(||V„v||| + h^^WQb^o - MIt) 



which gives rise to (|A.5[) after a summation over all T e T/j. □ 

Lemma A. 3. Let 1 < r < k and w (E [H''+^{n)]'^ and p G i?''(f2) and v e Vh- 
Assume that the finite element partition Th is shape regular. Then, the following 
estimates hold true 



(A.6) 
(A.7) 
(A.8) 



|s(Q?.w,v)| <C/in|w||.+i|||vI|l, 

i^w(v)i <c/iniwii.+iiiiv||i, 

|0p(v)| <C/inip||r|||v|||, 



where ^w(') o,nd £p{-) are two linear functionals on Vh given by 

(A.9) £w(v) = ^ (vo - Vb, Vw • n - Q,,(Vw) • n)^^,, 

TeTh 

(A.IO) 0p(v) = ^ (vo - vf,, {p - QhP)n)9T- 

Proof. Using the definition of Qi,, (jA.4[) and (jA.l[) . we have 



|s(Q/iW, v)| = 



dT 



Ten 



/ij.^((3b(Qow - w), QbVo-Vb)aT 



Ten 



dT 

Ten 

< I J2 (/^t'IIQow - will + ||V(gow - w)l||) 
\Ten J 

1/2 



1/2 



< C/i'-||w|L+i|||v|| 
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It follows from (|A.4|) and (|A.2|) that 



|C(v)| = 



^ (vo - Vb, Vw • n - Q;j(Vw) • n) 



dT 



< 



(vo - Qbvo, Vw • n - Q/i(Vw) • n) 



OT 



Ten 



^ (QfcVo - Vfc, Vw • n - Q/i(Vw) ■ n) 



OT 



Ten 



To estimate the first term on the righ-hand side of the above inequahty, we use (jA.4[) . 
(|A.2[) . (|A.5P and the inverse inequality to obtain 



^ (vo - QbVo, Vw • n - Qh(Vw) • n) 



dT 



Ten 



<C J2 hrW^yv ■ n - Q,,(Vw) ■ n||aT||Vvo| 



OT 



Ten 



\Ten J \Ten 

<c/iniw||^+iiiiv|||. 

Similarly, for the second term, we have 

— Vft, Vw • n — Qh(Vw) • n)^;j. 



1/2 



Vv( 



OWdT 



Ten 



<c(y1 /iT||Vwn-Q^(Vw)-n|||y^ ( hr^\\Qb^o-^b 
\Ten / \Ten 



1/2 



I ST 



<Ch'' 



|W||r+l|||V|| 



The estimate (IA.7P is verified by combining the above three estimates. 

The same technique for proving (jA.7l) can be applied to yield the following esti- 
mate. 



= vfc, ip - Qhp)n)9T 

Ten 

< c/inipiiriivii. 



This completes the proof of the lemma. □ 
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